﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.NumericalIntegration
{
    /// <summary>
    /// This class provides the numerical integration of integrals by using the Gauss-Kronod method.
    /// </summary>
    [Serializable]
    public class GaussKronrodIntegrator : AbstractGaussKronrodIntegrator
    {
        /// <summary>
        /// Initializes a new instance of the <see cref="GaussKronrodIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        /// <param name="integrationMethod">The integration method to use.</param>
        public GaussKronrodIntegrator(double a, double b, IRealFunction function,
                                      GaussKronrodIntegrationMethod integrationMethod)
            : base(a, b, function, integrationMethod)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="GaussKronrodIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        /// <param name="integrationMethod">The integration method to use.</param>
        public GaussKronrodIntegrator(int a, int b, SimplePolynomial function,
                                      GaussKronrodIntegrationMethod integrationMethod)
            : base(a, b, function, integrationMethod)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="GaussKronrodIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        /// <param name="integrationMethod">The integration method to use.</param>
        public GaussKronrodIntegrator(int a, int b, Polynomial function, GaussKronrodIntegrationMethod integrationMethod)
            : base(a, b, function, integrationMethod)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="GaussKronrodIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        /// <param name="integrationMethod">The integration method to use.</param>
        public GaussKronrodIntegrator(int a, int b, SimpleRationalPolynomial function,
                                      GaussKronrodIntegrationMethod integrationMethod)
            : base(a, b, function, integrationMethod)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="GaussKronrodIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        /// <param name="integrationMethod">The integration method to use.</param>
        public GaussKronrodIntegrator(int a, int b, HardRealFunction function,
                                      GaussKronrodIntegrationMethod integrationMethod)
            : base(a, b, function, integrationMethod)
        {
        }

        /// <summary>
        /// Integrates the integral by using the Gauss-Kronrod method.
        /// </summary>
        /// <param name="epsabs">The absolute requested accuracy of the result.</param>
        /// <param name="epsrel">The relative requested accuracy of the result.</param>
        /// <returns>The area of the specified integral.</returns>
        public double Integrate(double epsabs, double epsrel)
        {
            const int iterations = 500;
            const int limit = iterations - 1;
            int irule = (int) this.IntegrationMethod;
            double c = 0;
            double area1 = 0;
            double area2 = 0;
            double defab1 = 0;
            double defab2 = 0;
            double error1 = 0;
            double error2 = 0;
            double defabs = 0;
            double resabs = 0;
            double result = 0;
            int[] iord = new int[iterations];
            double[] alist = new double[iterations];
            double[] blist = new double[iterations];
            double[] rlist = new double[iterations];
            double[] elist = new double[iterations];

            int ier = 0;
            int neval = 0;
            int last = 0;
            double abserr = 0.0;
            alist[0] = this.A;
            blist[0] = this.B;
            rlist[0] = 0.0;
            elist[0] = 0.0;
            iord[0] = 0;

            if ((epsabs < 0.0) && (epsrel < 0.0))
            {
                ier = 6;
            }

            if (ier == 6)
            {
                this.absoluteAccuracy = epsabs;
                this.relativeAccuracy = epsrel;
                this.absoluteError = abserr;
                this.neededFunctionEvaluations = neval;

                return result;
            }

            int keyf = irule;

            if (irule <= 0)
            {
                keyf = 1;
            }

            if (irule >= 7)
            {
                keyf = 6;
            }

            c = keyf;
            neval = 0;

            switch (keyf)
            {
                case 1:
                    result = GaussKronrodRules.G_K15(this.Function, this.A, this.B, ref abserr, ref defabs, ref resabs);
                    break;
                case 2:
                    result = GaussKronrodRules.G_K21(this.Function, this.A, this.B, ref abserr, ref defabs, ref resabs);
                    break;
                case 3:
                    result = GaussKronrodRules.G_K31(this.Function, this.A, this.B, ref abserr, ref defabs, ref resabs);
                    break;
                case 4:
                    result = GaussKronrodRules.G_K41(this.Function, this.A, this.B, ref abserr, ref defabs, ref resabs);
                    break;
                case 5:
                    result = GaussKronrodRules.G_K51(this.Function, this.A, this.B, ref abserr, ref defabs, ref resabs);
                    break;
                case 6:
                    result = GaussKronrodRules.G_K61(this.Function, this.A, this.B, ref abserr, ref defabs, ref resabs);
                    break;
            }

            last = 0;
            rlist[0] = result;
            elist[0] = abserr;
            iord[0] = 0;

            double errbnd = (((epsabs) > (epsrel * Math.Abs(result))) ? (epsabs) : (epsrel * Math.Abs(result)));

            if ((abserr <= 50.0 * Double.Epsilon * defabs) && (abserr > errbnd))
            {
                ier = 2;
            }

            if ((ier != 0) || ((abserr <= errbnd) && (abserr != resabs)) || (abserr == 0.0))
            {
                if (keyf != 1)
                {
                    neval = (10 * keyf + 1) * (2 * (neval) + 1);
                }
                else
                {
                    neval = 30 * (neval) + 15;
                }

                this.absoluteAccuracy = epsabs;
                this.relativeAccuracy = epsrel;
                this.absoluteError = abserr;
                this.neededFunctionEvaluations = neval;

                return result;
            }

            double errmax = abserr;
            int maxerr = 0;
            double area = result;
            double errsum = abserr;
            int nrmax = 0;
            int iroff1 = 0;
            int iroff2 = 0;

            for (last = 1; last <= limit; (last)++)
            {
                double a1 = alist[maxerr];
                double b1 = 0.5 * (alist[maxerr] + blist[maxerr]);
                double a2 = b1;
                double b2 = blist[maxerr];

                switch (keyf)
                {
                    case 1:
                        area1 = GaussKronrodRules.G_K15(this.Function, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = GaussKronrodRules.G_K15(this.Function, a2, b2, ref error2, ref resabs, ref defab2);
                        break;

                    case 2:
                        area1 = GaussKronrodRules.G_K21(this.Function, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = GaussKronrodRules.G_K21(this.Function, a2, b2, ref error2, ref resabs, ref defab2);
                        break;

                    case 3:
                        area1 = GaussKronrodRules.G_K31(this.Function, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = GaussKronrodRules.G_K31(this.Function, a2, b2, ref error2, ref resabs, ref defab2);
                        break;

                    case 4:
                        area1 = GaussKronrodRules.G_K41(this.Function, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = GaussKronrodRules.G_K41(this.Function, a2, b2, ref error2, ref resabs, ref defab2);
                        break;

                    case 5:
                        area1 = GaussKronrodRules.G_K51(this.Function, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = GaussKronrodRules.G_K51(this.Function, a2, b2, ref error2, ref resabs, ref defab2);
                        break;

                    case 6:
                        area1 = GaussKronrodRules.G_K61(this.Function, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = GaussKronrodRules.G_K61(this.Function, a2, b2, ref error2, ref resabs, ref defab2);
                        break;
                }

                neval += 1;

                double area12 = area1 + area2;
                double erro12 = error1 + error2;

                errsum = errsum + erro12 - errmax;
                area = area + area12 - rlist[maxerr];

                if ((defab1 != error1) && (defab2 != error2))
                {
                    if ((Math.Abs(rlist[maxerr] - area12) <= 1.0e-5 * Math.Abs(area12)) && (erro12 >= .99 * errmax))
                    {
                        iroff1++;
                    }

                    if ((last > 9) && (erro12 > errmax))
                    {
                        iroff2++;
                    }
                }

                rlist[maxerr] = area1;
                rlist[last] = area2;
                errbnd = (((epsabs) > (epsrel * Math.Abs(area))) ? (epsabs) : (epsrel * Math.Abs(area)));

                if (errsum > errbnd)
                {
                    if ((iroff1 > 6) || (iroff2 > 20))
                    {
                        ier = 2;
                    }

                    if (last == limit)
                    {
                        ier = 1;
                    }

                    if ((((Math.Abs(a1)) > (Math.Abs(b2))) ? (Math.Abs(a1)) : (Math.Abs(b2))) <=
                        (1.0 + c * 1000.0 * Double.Epsilon) * (Math.Abs(a2) + 1.0e4 * Double.MinValue))
                    {
                        ier = 3;
                    }
                }

                if (error2 <= error1)
                {
                    alist[last] = a2;
                    blist[maxerr] = b1;
                    blist[last] = b2;
                    elist[maxerr] = error1;
                    elist[last] = error2;
                }
                else
                {
                    alist[maxerr] = a2;
                    alist[last] = a1;
                    blist[last] = b1;
                    rlist[maxerr] = area2;
                    rlist[last] = area1;
                    elist[maxerr] = error2;
                    elist[last] = error1;
                }

                DqSort(limit, last, ref maxerr, ref errmax, elist, iord, ref nrmax);

                if ((ier != 0) || (errsum <= errbnd))
                {
                    break;
                }
            }

            result = 0.0;

            for (int k = 0; k <= last; k++)
            {
                result += rlist[k];
            }

            abserr = errsum;

            if (keyf != 1)
            {
                neval = (10 * keyf + 1) * (2 * (neval) + 1);
            }
            else
            {
                neval = 30 * (neval) + 15;
            }

            this.absoluteAccuracy = epsabs;
            this.relativeAccuracy = epsrel;
            this.absoluteError = abserr;
            this.neededFunctionEvaluations = neval;
            this.precisionError = !(ier == 0);
            return result;
        }

        private static void DqSort(int limit, int last, ref int maxerr, ref double ermax, double[] elist, int[] iord,
                                   ref int nrmax)
        {
            int i;
            int j;
            int isucc;

            if (last > 1)
            {
                goto _10;
            }

            iord[0] = 0;
            iord[1] = 1;
            goto _90;

            _10:
            double errmax = elist[maxerr];

            if (nrmax == 0)
            {
                goto _30;
            }

            int ido = (nrmax) - 1;

            for (i = 0; i <= ido; i++)
            {
                isucc = iord[nrmax - 1];

                if (errmax <= elist[isucc])
                {
                    goto _30;
                }

                iord[nrmax] = isucc;
                (nrmax)--;
            }

            _30:
            int jupbn = last;

            if (last > (limit / 2 + 2))
            {
                jupbn = limit + 3 - last;
            }

            double errmin = elist[last];
            int jbnd = jupbn - 1;
            int ibeg = nrmax + 1;

            if (ibeg > jbnd)
            {
                goto _50;
            }

            for (i = ibeg; i <= jbnd; i++)
            {
                isucc = iord[i];

                if (errmax >= elist[isucc])
                {
                    goto _60;
                }

                iord[i - 1] = isucc;
            }

            _50:
            iord[jbnd] = maxerr;
            iord[jupbn] = last;
            maxerr = iord[nrmax];
            ermax = elist[maxerr];

            return;

            _60:
            iord[i - 1] = maxerr;

            int k = jbnd;

            for (j = i; j <= jbnd; j++)
            {
                isucc = iord[k];
                if (errmin < elist[isucc])
                {
                    iord[k + 1] = last;
                    maxerr = iord[nrmax];
                    ermax = elist[maxerr];

                    return;
                }

                iord[k + 1] = isucc;
                k--;
            }

            iord[i] = last;

            _90:
            maxerr = iord[nrmax];
            ermax = elist[maxerr];

            return;
        }
    }
}